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ABSTRACT 


A  new  method  of  array  modeling  which  will  be  used  to  predict 
the  performance  of  low  frequency  active  sonar  arrays  is  being 
investigated.  In  support  of  this  effort,  a  network  representation  of 
a  spherical  shell  piezoelectric  transducer  was  developed.  The 
transducer  was  modeled  using  the  finite-element  code  MARTS  AM, 
from  which  a  nodal  description  of  the  transducer  was  obtained.  A 
procedure  was  developed  to  reduce  and  transform  the  nodal 
description  of  the  transducer  into  a  spherical  harmonic  description. 
The  spherical  harmonic  description  of  the  transducer  was  computed 
at  two  frequencies,  112.5  Hz  and  1125.3  Hz,  corresponding  to  values 
of  ka  of  0.1  and  1.0,  respectively,  where  a  is  the  radius  of  the 
sphere. 
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I.  INTRODUCTION 


The  objectives  of  the  research  detailed  in  this  thesis  are:  1)  to 
produce  a  nodal  description  of  a  spherical  shell,  piezoelectric 
transducer  using  finite-element  modeling,  and  2)  to  develop  a 
procedure  to  transform  the  nodal  description  of  the  transducer  into 
a  spherical  harmonic  description. 

This  report  is  organized  in  three  main  sections.  First,  a  new 
method  to  predict  the  performance  of  low  frequency  active  arrays, 
which  provided  the  motivation  for  this  research,  is  described.  Next, 
the  finite-element  modeling  and  description  in  terms  of  spherical 
harmonics  of  a  particular  spherical  shell  transducer  are  detailed. 
Numerical  results  are  included  for  the  transducer  operating  at  112.5 
Hz  and  1125.3  Hz,  corresponding  to  values  of  ka  of  0.1  and  1.0, 
respectively,  where  a  is  the  radius  of  the  sphere.  Lastly,  the 
procedure  to  reduce  and  transform  the  nodal  description  of  the 
transducer  to  a  spherical  harmonic  description  is  described. 
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II.  MODELING  OF  LOW  FREQUENCY  ACTIVE  ARRAYS 


A.  BACKGROUND 

The  trend  toward  operation  of  sonar  surveillance  systems  at 
lower  frequencies  has  necessitated  larger,  more  dense  arrays  of 
transducers.  This  transition  has  prompted  investigation  of  a  new, 
potentially  more  efficient,  method  to  predict  array  performance. 
This  new  method  of  array  modeling  combines  a  finite*element 
representation  for  each  transducer  with  a  mathematical 
representation  of  the  acoustic  field  which  is  equivalent  to  the 
T-matrix  method  which  has  been  applied  to  elastic  scattering 
problems  [Ref.  1].  This  approach  is  explained  in  greater  detail  in 
Appendix  A.  It  has  the  feature  that  it  accounts  for  the  effects  of 
'multipie-scattering',  which  can  be  quite  significant  in  dense  arrays 
(the  distortion  of  the  near-field  radiation  pattern  of  a  transducer 
due  to  a  nearby  transducer,  for  example,  is  a  manifestation  of 
multiple-scattering).  Applying  this  method  to  an  array  of 
transducers  will  ultimately  yield  a  matrix  which  relates  the 
outgoing  acoustic  waves  to  the  input  electrical  excitation  of  each 
transducer. 

Figure  1  is  a  schematic  diagram  of  a  portion  of  an  array.  In 
this  method,  the  modeling  problem  is  divided  into  two  regions.  The 
term  structure  refers  to  the  transducer  plus  some  surrounding  fluid 
out  to  an  imaginary,  spherical  surface,  S.  The  second  region  is  the 
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fluid  outside  of  S.  The  array  may  be  composed  of  any  nurrSer  and/or 
variety  of  tranducers.  Additionally,  it  may  have  any  geometry. 

The  goal  is  to  develop  an  analytical  representation  of  the 
structure  from  a  finite-element  model  and  to  combine  it  with  an 
analytical  representation  of  the  acoustic  field  in  the  fluid  to  predict 
the  performance  of  an  array  of  interacting  transducers.  By 
separately  modeling  the  structure  and  the  fluid,  this  approach  will 
simplify  the  complex  task  of  array  modeling.  Changing  the  geometry 
of  an  array  will  require  recomputing  only  the  fluid  model.  Likewise, 
maintaining  the  geometry  and  changing  the  type  of  transducer  will 
require  a  new  representation  of  only  the  sructure. 

B.  FLUID  REPRESENTATION 

The  acoustic  pressure  field  in  the  unbounded  fluid  is 
represented  as  an  eigenfunction  expansion  in  terms  of  spherical 
waves.  It  should  be  noted  that  all  quantities  are  assumed  to  vary 
harmonically  with  time.  This  representation  is  convenient  due  to 
the  spherical  shape  of  the  structures.  Figure  2  is  the  radiation 
pattern  which  is  produced  by  three  pulsating  spherical  elastic  shells 
of  radius  a,  with  one  half  wavelength  separation  and  ka  equal  to  1.0. 
Each  shell  is  driven  with  a  uniform  pressure  amplitude  on  the  inside 
of  the  surface  of  the  strength  indicated  by  the  number  inside.  The 
curves  shown  represent  the  radiation  pattern  when  the  number  of 
spherical  harmonics  retained  in  the  expansion  is  varied.  The  curves 
which  include  the  first  6,  9,  and  11  harmonics  are  indistinguishable: 
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Figure  2:  Radiation  Pattern  of  Three  Spherical  Shells. 
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therefore,  the  acoustic  field  may  be  accurately  described  by 
retaining  only  about  6  spherical  harmonics.  [Ref.2] 

C.  STRUCTURE  REPRESENTATION 

The  structure  must  be  spherical  in  shape;  it  may  be  a 
transducer  with  some  surrounding  fluid  or  simply  a  spherical 
transducer.  The  structure  is  first  modeled  using  a  finite-element 
code,  which  produces  a  nodal  description  of  the  structure.  The  nodal 
description  is  then  reduced  and  transformed  into  a  spherical 
harmonic  description  of  the  structure.  This  description  can  then  be 
combined  with  the  spherical  harmonic  description  of  the  acoustic 
field  on  the  surface  of  each  sphere  to  predict  the  behavior  of  the 
array. 
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III.  MODELING  OF  SPHERICAL  SHELL  TRANSDUCER 


A.  FINITE  ELEMENT  MODEL 

A  spherical  shell  piezoelectric  transducer  was  modeled  using 
the  finite  element  code  MARTSAM.  Figure  3  is  the  finite  element 
mesh  used  to  model  the  transducer.  The  transducer  is  radially 
polarized  Navy  Type  1  ceramic  with  an  outer  radius  of  7.620 
centimeters  and  thickness  of  0.762  centimeters.  It  is  symmetric 
about  the  y-axis.  Based  on  this  symmetry,  it  was  only  necessary  to 
create  a  two-dimensional  model  of  one  half  of  the  transducer.  The 
structure  was  divided  into  60  six-noded  triangular  elements.  This 
partitioning  scheme  results  in  a  mesh  with  183  nodes.  There  are 
two  mechanical  degrees  of  freedom  associated  with  each  node 
except  the  six  nodes  along  the  y-axis.  Boundary  conditions  were 
applied  to  set  the  motion  of  these  nodes  in  the  x-direction  equal  to 
zero.  Consequently,  this  structure  possesses  a  total  of  360 
mechanical  degrees  of  freedom.  Additional  boundary  conditions 
included  grounding  the  electrode  on  the  outer  surface  of  the  sphere 
while  maintaining  a  constant  potential  on  the  inner  electrode.  The 
structure,  therefore,  has  one  electrical  degree  of  freedom. 

B.  MARTSAM  OUTPUT 

A  modal  analysis  of  the  structure  was  performed  by  Dr.  Michele 
McCollum  of  the  Naval  Research  Laboratory,  Orlando  using  the  finite- 
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Figure  3:  Finite  Element 
Transducer. 


of  Spherical  Shell 


element  code  MARTSAM.  For  this  particular  application,  the  code 
was  used  to  generate  (Kuu)>  (M^,  and  Kuv  which  can  be  substituted 

into  the  following  set  of  equations: 

where 

•  ii  and  F  are  vectors  that  contain  the  nodal  values  of 
the  displacement  field  and  the  applied  forces, 

•  V  is  the  applied  electrical  potential, 

•  Q  is  the  electrical  charge  on  the  structure, 

•  (Kuu)  is  the  matrix  which  describes  the  effect  that  the 

stiffness  at  each  node  has  on  alt  the  nodes, 

.(M)  is  the  matrix,  which  describes  the  effect  that  the 

mass  at  each  node  has  on  all  the  nodes, 

•  j<uv  is  the  vector  that  contains  the  coupling  coefficients 
that  relate  the  mechanical  and  electrical  degrees  of 
freedom, 

•  Kvv  is  the  capacitance  of  the  transducer  for  zero 
displacement  everywhere, 

•  CO  is  the  angular  frequency, 

•  T  means  transpose. 
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to  completely  describe  the  spherical  shell  transducer  in  terms  of  its 
nodes  at  a  specified  frequency.  The  matrix  of  the  left  side  of 
equation  (3.1)  has  dimension  361  by  361  for  this  structure. 

By  means  of  matrix  algebra,  which  is  outlined  in  Chapter  IV, 
this  nodai  description  of  the  structure  was  reduced.  The  mechanical 
degrees  of  freedom  associated  with  the  nodes  internal  to  the 
structure  were  eliminated  as  the  external  force  on  those  nodes  is 
zero.  The  displacements  at  the  surface  nodes  were  transformed 
from  the  Cartesian  coordinate  system  to  the  poiar  coordinate 
system.  The  degrees  of  freedom  associated  with  the  parallel 
displacements  were  also  eliminated  as  external  fluid  forces  are 
applied  normal  to  the  surface.  The  dimensions  of  the  reduced  matrix 
are  62  by  62,  61  surface  normal  displacement  degrees  of  freedom 
and  1  electrical  degree  of  freedom. 

C.  SPHERICAL  HARMONIC  DESCRIPTION  OF  STRUCTURE 

Following  the  procedure  detailed  in  Chapter  IV,  the  reduced 
nodal  description  of  the  transducer  was  transformed  into  a  spherical 
harmonic  description.  Since  the  finite-element  model  is  restricted 
to  axisymmetric  solutions  for  field  quantities,  the  Legendre 
polynomials  were  used  for  the  spherical  harmonic  functions.  Based 
on  the  results  obtained  by  Canright  and  Scandrett  [Ref.  2],  the  first  7 
Legendre  polynomials  were  used.  The  spherical  harmonic  description 
of  the  structure  is  an  8  by  8  matrix.  Figures  4  and  5  are  the 
electrodynamic  matrices  obtained  for  the  frequencies  112.5  and 
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1125.3  Hz,  respectively,  corresponding  to  values  of  ka  of  0.1  and  1.0, 

respectively,  where  a  is  the  radius  of  the  structure.  These 
matrices,  are  used  in  the  following  set  of  equations: 

(■"KD-pr).  M 

where  and  are  vectors  that  contain  the  values  of  the 
displacement  field  and  stress  field  in  terms  of  Legendre  polynomials 
(n-0,1,2...6),  to  describe  the  transducer.  Most  of  the  elements  in 
these  matrices  are  on  the  order  10®  -  lO"*®.  These  numbers 
represent  mechanical  interactions  while  the  elements  in  the  last 
row  and  column  (on  the  order  of  10'®  -  10)  represent  the  mechanical 
and  electrical  coupling.  The  element  in  the  last  row  and  column  is 
the  (surface  only)  blocked  capacitance. 
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for  a  Spherical  Shell  Transducer,  f  =  112,5. 
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for  a  Spharical  Shell  Transducer,  f  =  1125.3. 


IV.  TRANSFORMATION  PROCEDURE 


A.  REDUCED  NODAL  DESCRIPTION 

Matrix  equation  (3.1),  for  simplicity,  can  be  represented  by 

(pnod)  i^nod  *  pnod,  (4,1) 


where 

.  (D~a)  is  the  square  dynamical  matrix. 

•  y^nod  is  the  vector  containing  the  displacement  field 
and  the  electrical  potential,  and 

•  F^od  is  the  vector  containing  the  applied  force  field 
and  the  electrical  charge. 

The  dimensions  of  the  dynamical  matrix  are  quite  large  for  even  the 
simplest  transducer.  In  modeling  the  fluid-loaded  structure, 
however,  the  size  of  this  description  can  be  greatly  reduced  by 
applying  matrix  algebra  as  the  following  procedure  describes.  In 
performing  this  reduction,  the  accuracy  of  the  description  is  not 
diminished  since  the  fluid  forces  act  only  at  the  surface  bounding 
the  structure. 

1.  Forces  and  Internal  Nodes 

£nod  is  initially  a  column  vector  whose  entries  are  the  force 
applied  to  each  node  listed  by  ascending  node  number  and  direction 
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followed  by  the  electrical  charge.  However,  the  applied  force  is 
zero  at  ail  the  nodes  that  are  internal  to  the  structure  as  those 
nodes  are  inaccessible  to  the  surrounding  fluid.  is  rearranged  so 
that  all  the  zero  entries  come  first.  When  ,  the  vector  on  the 
right  hand  side  of  the  equation,  is  reordered,  the  rows  of  the 
dynamical  matrix  on  the  left  side  of  the  equation  must  be  rearranged 
accordingly.  must  also  be  reordered  in  the  same  manner  as 

which  necessitates  a  corresponding  rearranging  of  the  columns  of 
the  dynamical  matrix. 

2.  Transformation  of  Coordinate  System 

MARTSAM  describes  the  transducer  in  terms  of  Cartesian 
coordinates.  A  transformation  matrix,  can  be  defined  for  each 

surface  node  such  that  it  performs  the  following  operations; 


where,  respectively, 

•  Ux  and  Uy  are  the  displacements  in  the  x-  and  y-directions, 

•  £x  and  Ey  are  the  applied  forces  in  the  x-  and  y-directions, 

•  Uperp  and  Upar  are  the  displacements  in  the  directions 
perpendicular  and  parallel  to  the  surface,  and 
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•  Eperp  and  Epar  are  the  applied  forces  in  the  directions 
perpendicular  and  parallel  to  the  surface. 

For  the  mesh  shown  in  Figure  3,  for  a  single  element, 

,  ,  f  sine  -cose  r  X  -Y  ^ 

n-i  cose  Sine  Y  X  )•  (“-a) 

where 

•  6  is  the  angle  each  surface  node  makes  with  the  positive 
y-axis, 

•  X  and  Y  are  the  Cartesian  coordinates  of  each  node,  and 

•  R  is  the  radially  distance  of  each  node  from  the  center  of 
the  sphere. 

Substituting  equations  (4.2)  into  equation  (4.1)  yields 

(Onod)  ^T)  =  (T)  EPo'ar,  (4.4) 

where  the  vectors  U  and  F  are  now  expressed  in  terms  of  polar 

coordinates.  Multiplying  each  side  of  equation  (4.4)  by  the  inverse  of 
(T)  gives 


(T)-1  (Drod)  (T)  iiPolar  s£polar^ 


(4.5) 
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where  (T)‘^  (*^"**)  (T)  represents  the  dynamical  matrix  in  terms  of 

nodes  and  the  polar  coordinate  system. 

3.  Forces  and  Parallel  Direction 

The  elements  of  £Poiar  are  arranged  as  follows:  the  forces 
applied  to  the  internal  nodes,  the  forces  applied  in  the  directions 
perpendicular  to  the  transducer  at  the  surface  nodes  and  parallel  to 
the  transducer  at  the  surface  nodes,  and  the  electrical  charge. 
However,  since  only  the  forces  applied  normal  to  the  transducer  are 
non-zero,  Fpar  is  the  null  vector.  Using  the  procedure  described  in 
section  1 ,  Epolar  and  UPolar  must  be  reordered  such  that  Epar  and  iipar 
precede  Fperp  and  Uperp.  The  rows  and  columns  of  the  dynamical 
matrix  must  be  rearranged  accordingly. 

4.  Reduction  of  Dynamical  Matrix 

The  procedures  in  sections  1,  2  and  3  have  arranged  the 
dynamical  matrix  into  the  proper  form  to  be  reduced  using  matrix 
algebra.  The  matrix  equation  is  now  in  the  form 

(111}  !1}r"}  )(uSJ=(epIp).  <-) 


where 

•  (U_),(UR),(LL),(LR)  are  subdivisions  of  the 
dynamical  matrix, 

•  Uo  is  the  displacements  at  the  internal  nodes  and  in  the 
direction  parallel  to  the  surface  nodes, 
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*  iiperp  is  the  normal  displacements  and  the  electrical 
potential, 

*  0.  is  the  null  vector,  and 

*  Eperp  are  the  forces  applied  normal  to  the  surface  and  the 
electrical  charge. 


The  corresponding  simultaneous  equations  are 

(UL)  iio  +  (i^)  iiperp  “  fil .  and  (4.7) 

(l-i-)  iio  +  iiperp  *  Eperp  ■  (4.8) 

Solving  equation  (4.7)  for  Uo  yields 

iio  =  -  (UL)-1  (UR)Uperp.  (4.9) 

Substituting  equation  (4.9)  into  equation  (4.8)  gives 

((m)-(LL)(UL)-l(m))g^,,p.  Eperp,  (4.10) 


where  ((^)"  (^)  (^)’ '  {*^) )  '®  reduced  dynamical  matrix 

describing  nodal  interactions.  This  matrix  is  square,  and  its 
dimension  is  the  number  of  surface  nodes  plus  the  electrical  degree 
of  freedom.  The  computer  code,  which  computes  the  reduced 
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dynamical  matrix  for  the  spherical  shell  transducer  described  in 
Chapter  III,  is  listed  in  Appendix  B  [Ref.  3:p.  38]. 

B.  TRANSFORMATION  TO  SPHERICAL  HARMONIC  DESCRIPTION 
The  reduced  nodal  description  of  the  transducer  given  by 
equation  (4.10)  can  be  written  as 


'  (DSS") 

> 

Knod 

^piod.. 

KnodT 

Cffi  ^ 

)l  V  J  = 

‘1  °J 

where 

•  (Duu*^)  is  the  reduced  (Kuu)  -  matrix, 

•  l^nod  is  the  reduced  vector  containing  the  coupling 
coefficients  that  relate  the  mechanical  and  electrical 
degrees  of  freedom, 

•  Ceb  is  the  (surface  only)  blocked  capacitance, 

•  y_nod  is  the  vector  containing  normal  displacements, 

•  V  is  the  electrical  potential, 

•  £nod  is  the  vector  containing  the  forces  applied  normal  to 
the  surface, 

•  Q  is  the  electrical  charge,  and 

•  T  means  transpose. 

A  spherical  harmonic  description  of  the  transducer  can  be  obtained, 
in  the  form 


1  9 


(DsphJ  i£sph  =  Qsph  , 


(4.12) 


by  performing  the  following  tranformations. 

1.  Transformation  of  to  iL^ph 

The  surface  normal  displacement  distribution  can  be  written 
as 


N 

U(e)  =  SuR°‘‘  Nn(e)  ,  (4.13) 

n 

where  =  U(0  =  0n)  ,  and  N(0)  is  the  finite-element  interpolation 
function  where  Nn(0)  =  5(0  =  0n)  [Ref.  4].  Since  the  finite-element 
model  created  to  represent  the  trarscucer  in  this  investigation  is 
restricted  to  axisymmetric  solutions  for  all  field  quantities,  the 
Legendre  polynomials  were  chosen  for  the  set  of  harmonic  functions. 
The  surface  normal  displacement  distribution  can  then  be  expressed 
approximately  by  a  truncated  series  as 

M 

U(0)  =  Pm(cos0)  .  (4.14) 

m 


n 

where  -  (1/Cm)  Ju(0)  Pm(cos0)  sin0  d0  and  Cm  =  J(Pm)2  sin0  dS, 

0 

the  usual  normalization  constant  associated  with  Pm.  For  the 
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Legendre  polynomials,  Cm  «  (21-1-1  )/2  where  I  is  the  order.  From 
equations  (4.13)  and  (4.14),  it  follows  that 

M 

=XUm^^  Pm(cos0=cos0n)  -  (4.15) 

m 

Equation  (4.15)  expressed  in  matrix  form  is 

iinod«  (p)U8ph,  (4  16) 

where  the  matrix  elements  Pnm  =  Pm(cos0=cos0n)  . 

2.  Transformation  of  f.nod  to 

The  self-consistent  force  at  node  n  due  to  a  distributed 
stress  o  is  given  by 


Fn'»«  =  JoNndS, 


(4.17) 


where  the  integration  is  performed  over  the  surface  of  the 
structure.  The  distributed  stress  can  be  expanded  similarly  to  the 
displacement  as 


M 

0(0)  =  Pm(cos0) .  (4.18) 

m 
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The  dynamical  equation  describing  the  transducer  in  terms  of  s^Ph 
follows  from  the  equation  for  the  coefficients  of  the  expansion  of 
the  stress  field  in  the  Legendre  polynomials: 

*  (1/Cm)  Jo  Pm  sine  dS.  (4.19) 

The  expansion. 

N 

Pm(e)  =  Ip  nm  Nn  ,  (4.20) 

n 

where  Pnm  Pm(9n).  can  be  made.  Substitution  of  equation  (4.20) 
into  equation  (4.19)  gives 

h  N 

-(1/Cm)I  (J<jN„dS)  Pmn.  (4.21) 

n 

Using  equation  (4.17),  equation  (4.21)  becomes 

N 

=  (1/Cm)  Pmn .  (4-22) 

n 

In  matrix  form,  equation  (4.22)  can  be  written 

ffsph  =  (C)-1  (  P)1’  pnod  ,  (4.23) 
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where 


•  (C)  =  diag  (Cm) , 

•  (  P)nm  =  Pm(cosen)  •  and 

•  FR°‘‘  =  JoNndS. 


Equation  (4.23)  is  the  desired  expression  relating  ffSPh  to  Enod. 
Matrix  equation  (4.1 1 )  can  be  expressed  as 


(dGS**)  i^nod  +  i^nod  v  =  £nod  ,  and 

(4.24) 

j^nod  T  ^nod  +  Ceb  V  =  Q  . 

(4.25) 

Substitution  of  equation  (4.24)  into  equation  (4.23)  gives 

ffsph  =  (C)  - 1  (  P)  ■*■  ( +  (^) ‘  ^  (  P) ^  •  ("^-26) 

Therefore,  by  substituting  equation  (4.16),  equation  (4.26)  reduces 
to 

(d56^)  iisph  +  j^sph  V  «  ffsph  ,  (4.27) 

where  (d55^)  =  (C)-i  (P)T  (Duu^)  (P)  andK®P^  =  (C)-^  (P)Tlcnod. 

It  follows  from  equations  (4.16)  and  (4.25)  that 

l^nod  T  ^p)  j^sph  +  Ceb  V  =  Q  .  (4.28) 
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Hence  the  reduced  dynamical  matrix  equation  in  terms  of 
Legendre  polynomials  coefficient  degrees  of  freedom  is 


(C)-1  (  P)T  K.nod  "j^sphx 

linodT(p)  "  '  Jl  V  )-[  Q  J- 


(4.29) 


The  computer  code,  which  computes  the  transformation  from  the 
reduced  nodal  dynamical  matrix  to  the  spherical  harmonic  matrix,  is 
listed  in  Appendix  C  [Ref.  3:p.  38]. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  represents  the  first  step  in  an  investigation  which, 
when  successfully  completed,  will  significantly  improve  the  Navy's 
ability  to  predict  the  performance  of  dense  sonar  arrays.  A  network 
(matrix)  representation  of  a  spherical  shell  transducer  has  been 
calculated  in  terms  of  spherical  harmonics.  The  procedure  for 
reducing  and  transforming  a  nodal  description  of  a  transducer  into  a 
spherical  harmonic  description  has  been  documented. 

For  this  thesis,  a  two-dimensional  model  of  a  spherical  shell 
transducer  was  created  and  a  procedure  was  developed  to  reduced 
and  transform  its  nodal  description  into  a  spherical  harmonic  one,  it 
is  recommended  that  a  three-dimensional  model  be  created  and  a 
similar  transformation  procedure  be  developed.  Since  it  is  desirable 
that  (o  remain  a  free  parameter  when  describing  a  transducer,  it  is 
suggested  that  the  procedure  detailed  in  Chapter  IV  be  rederived 
following  Benthien's  method  to  obtain  the  dynamical  matrix  as  a 
function  of  co  [Ref.  5].  It  is  also  recommended  that  experimental  data 
be  obtained  to  verify  the  matrix  representation  developed  in  this 
thesis. 
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APPENDIX  A: 

EXCERPTS  FROM  PROPOSAL  FOR  LFA  ARRAY  MODELING 


APPLICATION  OF  THE  T-MATRIX  METHOD  TO  LOW-FREQUENCY 
ACTIVE  ARRAY  PERFORMANCE  PREDICTION 

TECHNICAL  PROPOSAL 
A.  Introduction  and  Related  Woik 

The  direction  of  active  sonar  surveillance  systems  is  toward  lower  frequencies,  requiring 
arrays  of  large,  high  power  transducers.  The  successful  design  and  operation  of  such  arrays 
requires  the  ability  to  reliably  predict  their  performance. 

Array  performance  prediction  is  a  coupled  structure-fluid  problem  [1].  For  die  analysis  of 
complex  structures,  such  as  a  sonar  transducer,  the  finite-element  method  (FEM)  is  preeminent 
[2].  Many  FEM  conqiuter  codes  exist,  some  of  which  include  piezoelectric  elements  for  the 
analysis  of  piezoelectric  transducers,  e.g.  MARTSAM,  ATILA.  For  application  to  coupled 
structure-fluid  problems,  the  major  drawback  of  the  finite-element  method  is  the  difficulty  of 
modeling  the  irifinite  fluid.  So-called  boundary  elements,  derived  from  a  Helmholtz  Integral 
formulation  of  the  exterior  radiation  problem,  have  been  introduced  into  the  FEM  to  terminate 
the  fluid  mesh  [3,4,5].  These  elements  typically  inpose  an  outgoing  radiation  impedance 
conditiem  on  the  adjacent  fluid  element.  The  frequency  dependence  of  the  radiation  impedance 
is  commonly  approximated  by  an  interpolation  between  the  low-  and  high-frequency 
asyrrptotic  limits  [S].  The  a{plication  of  a  combined  FEM-boundary  integral  method  to  array 
performaiKe  prediction  has  not  a{peared  in  the  literature.  A  full  Hdmholtz  Integral  formulation 
can  be  rpplied  at  the  bounding  surface  [2,3],  but  Uiis  technique  is  limited  to  small  problems  by 
the  number  of  degrees  of  freedom,  and  so  is  not  practical  to  model  an  entire  array.  HeiKe  there 
is  a  need  for  economical  methods  for  array  performance  prediction.  This  proposal  addresses 
that  need. 

The  method  we  prepose  is  formulated  in  terms  of  coupled  networks,  me  for  each  radiating 
element  and  orK  for  die  acoustic  field.  It  is  equivalent  to  the  T-matrix  m^hod  which  has  been 
applied  to  scattering  proUems  [7,8,9,10],  in  that  ultimately  a  matrix  can  be  derived  which 
relates  the  outgoing  acoustic  wave  to  die  iiput  dectrical  excitation  of  each  transducer  (the 
equivalent  reciprocal  problem  is  to  relate  the.electrical  output  of  each  transducer  to  an  incident 
wave).  One  could  in  fact  consider  the  proposed  research  the  application  of  the  T-matrix 
method  to  sonar  array  performance  prediaion.  As  in  the  T-matrix  method,  we  adopt  the  idea 
of  representing  the  acoustic  field  as  an  eigenfimetion  expansion  and  solving  for  the  coefficients 
self-consistendy.  No  restriction  is  placed  on  die  arrangement  of  the  radiators,  and  so  multiple¬ 
scattering  of  all  orders  is  rigorously  included.  In  addition,  we  plan  to  explicidy  take  into 
account  the  dynamical  prrperties  of  the  radiators.  These  are  to  be  found  for  a  tea’  transducer 
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using  a  finite-element  computer  code  [2,13]. 
B.  Proposed  Effort 


Goal 

The  goal  of  the  proposed  research  is  to  produce  the  most  economical  yet  complete 
descr^tim  possible  of  sonar  array  performance,  witfi  q>ecific  application  to  low  fiequotQr 
active  arrays.  The  idea  is  to  combine  an  analytical  rq>resentation  of  a  single  element  (derived 
fitvn  a  finite-element  model  for  a  real  transducer)  wiA  an  analytical  representation  of  the 
acoustic  fidd  to  predict  the  perfoimance  of  an  array  of  interacting  transducers.  Especially  for 
the  usual  case  of  an  array  of  identical  elements,  since  the  dynamic  bdiavior  of  only  a  single 
element  need  be  computed,  and  since  this  computation  need  be  d<xie  (»ily  once,  regardless  of 
array  geometry,  a  variety  of  operating  frequencies  and  array  geometries  may  be  analyzed  most 
economically. 

Method 

A  schematic  diagram  of  a  portion  of  an  array  of  sonar  transducers  is  shown  below. 


Figure  1 .  Schematic  of  a  portion  of  an  array  of  sonar  transducers. 

The  transducers  need  not  be  identical,  nor  in  any  particular  arrangement.  There  may  or  may 
not  be  an  acoustic  wave  incident  from  a  source  external  to  the  array.  The  closed  surface  labeled 
S  indicates  the  boundary  between  "structure"  and  "fluid".  5  may  coincide  with  the  physical 
boundary  of  a  transducer,  or  it  may  include  some  surrounding  fluid. 
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There  are  several  reasons  to  choose  S  to  lie  within  the  fluid  surrounding  a  transducer.  The 
most  important  is  that,  if  the  "structure"-"fluid”  bouixlary  coincides  with  a  constant-coordinate 
surface  of  a  coordinate  system  in  which  the  wave  equation  is  separable,  dien  the  eigenfiinctims 
of  the  wave  equation  in  that  coordinate  system  are  a  particularly  convenient  choice  as  the  basis 
for  expanding  both  the  surface  normal  velocity  on  S  and  die  acoustic  field  within  die  "fluid". 
This  makes  it  much  easier,  for  exan^le,  to  compute  the  self-  and  mutual  radiation  inqiedances. 
Also,  by  iiKluding  some  surrounding  fluid  as  part  of  the  "stmcture",  the  burden  of  carrying 
many  higher-order  (higher  spatial  fiequency)  terms  in  an  eigenfimction  expansion  of  the 
acoustic  field,  which  might  be  required  to  describe  the  local  flow  field  near  a  real  transducer, 
would  be  relieved.  These  degrees  of  freedom  would  be  considered  as  internal  to  die  structure, 
with  the  result  that  far  fewer  degrees  of  fieedom  would  be  required  to  describe  the  acoustic 
field.  This  is  important  because  only  the  descri|Xion  of  the  acoustic  field  within  the  "fluid" 
needs  to  be  recalculated  for  a  change  in  array  geometry. 

Conversely,  by  including  some  surrounding  fluid,  the  description  of  the  acoustic  field 
would  be  unchanged  by  a  change  in  the  physical  stmcture,  for  example,  by  exchanging  one 
type  of  transducer  for  another  (e.g.  flextensional  for  hydroacoustic).  The  "stmcture"  could 
even  be  composed  of  more  than  one  physical  transducer.  This  would  be  useful,  for  example, 
to  analyze  the  performance  of  an  array  composed  of  groups  of  closely-spaced  transducers. 

Including  some  surrounding  fluid  as  part  of  the  "stmcture"  raises  the  possibility  of 
introducing  spurious  stmctural  resonances,  which  may  cause  problems  for  a  numerical 
computation  [10].  The  consequences  of  this  will  have  to  be  investigated. 

It  is  convenient  to  discuss  the  array  performance  problem  further  in  terms  of  a  coupled 
network  representation.  The  diagram  below  depicts  the  coupling  between  "stmcture"  and 
"fluid"  for  one  surface  normal  velocity  degree-of-fireedom.  For  simplicity,  the  transducer  is 
considered  to  be  reciprocal. 
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Surface  normal  velocity 
degree-of-freedom  n 

Transducer  i  ,,  Acoustic  field 

“in 


Figure  2.  Coupled  reciprocal  netivoik  {presentation  of  structure-fluid 
interaction  for  one  surface  normal  velocity  degree-of-freedom. 


The  network  equations  representing  transducer  i  and  the  acoustic  field  are 


Cj  = 

^E,ih 

-  X/ik“ik» 
k 

*  X2ni,ink**ik* 
k 

(transducer  i) 

/in  “ 

^inh 

/in  ~ 

DuAiPf+ 

X  ^ranjm^jm* 

(acoustic  field) 

j,m 


The  subscripts  i  and  j  refer  to  particular  transducers,  and  the  subscripts  k,  n  and  m  refer  to 
particular  surface  normal  velocity  degrees-of-ficedom  (DOF)  on  a  "structure" -"fluid" 
boundary,  assumed  to  be  a  constant-coordinate  surface  of  a  coordinate  system  in  which  the 
wave  equation  is  separable.  The  meaning  of  the  remaining  symbols  is 

=  the  voltage  across  the  electrical  terminals  of  transducer  i, 
ij  =  the  current  through  the  electrical  terminals  of  transducer  i, 

/in  =  the  modal  force  amplitude  of  transducer  i,  DOF  n, 

Ujn  =  the  modal  surface  normal  velocity  amplitude  of  transducer  i,  DOF  n, 

Z£  i  =  the  blocked  electrical  impedance  of  transducer  i. 
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^in,inn  “ 


T  = 
^  in 

Z  ■  = 


Din  = 


the  open-circuit  self  mechanical  inqiedance  of  transducer  i,  DOF  n, 

the  open-circuit  mutual  mechanical  in^iedance  between  transducer  i,  DOF  n,  and 

DOFm 

the  transduction  coefficient  of  transducer  i,  DOF  n, 
the  self  radiation  inqiedance  of  transducer  i,  DOF  n, 

the  mutual  radiation  in^)edance  between  transducer  i,  DOF  n,  and  transduco:  j, 
DOFm, 

the  incident  £tee-field  pressure,  if  any  (usually  assumed  to  be  a  plane  wave), 
the  surface  area  of  the  "structure"-"fluid"  boundary  of  transducer  i, 
the  diffraction  constant  of  transducer  i,  DOF  n  (the  diffraction  constant  equals  the 
ratio  of  the  blocked  modal  force  an^litude  per  unit  area  to  the  incident  finee-field 
pressure). 


The  acoustic  field  is  to  be  represented  as  the  superposition  of  the  acoustic  field  due  to  each 
radiator  and  an  optional  incidoit  field.  The  surface  normal  velocity  and  the  acoustic  field  of 
each  radiating  element  are  to  be  represented  as  expansions  in  the  free-space  eigoifunctions  of 
the  wave  equation.  ^  n  ddition  theorem  [8,1 1,12]  is  to  be  used  to  express  the  field  due  to  one 
radiator  at  the  surface  of  another.  The  ^in  ^  to  be  found  in  terms  of  the 

eigenfunction  expansions  by  applying  tl^  boundary  condition  at  the  surface  of  each  radiator 
that  Wjm=^j<^nn  “j,„=0,  respectively.  No  restriction  is  placed  on  the  arrangement  of 
the  radiators,  and  so  multiple-scattering  of  all  orders  is  rigorously  included.  A  self-consistent 
solution  for  the  coefficients  of  the  eigenfunction  expansions  for  each  radiator  is  to  be  found  by 
applying  an  impedance-matching  boimdary  condition  at  the  surface  of  each  radiator.  The 
mechanical  impedance  of  each  radiator  is  to  be  represented  in  terms  of  the  modal  parameters  of 
its  in-vacuo  eigenfunctions.  These  are  to  be  found  for  a  real  transducer  using  a  finite-element 
computer  code  [2,13]. 
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APPENDIX  B 


c  PROGRAM  REDDYNMTRX:  computes  the  reduced  nodal  dynamical  matrix 

c  for  LT  McLean’s  radially  polarized  piezoelectric  spherical  shell 

real  omega ,omega2, pi , f req, r r ,pint( 361 , 361 ) 
real  dyn( 361 , 361 ) , kuu( 360 , 360 ) ,mas( 360, 360) 
real  r ( 361 , 361 ) ,g( 361 , 361 ) 

real  a( 299 , 299 ) .b( 299,62 ). c( 62, 299) ,d( 62,62) 
real  nn(62,62) ,sum 

real  p(62, 299 ) ,pp( 62, 62 ) ,prod( 361, 361 ) 
real  s( 361 , 361 ) ,mat( 361 , 361 ) ,y( 299,299) , indx( 299 ) 
real  t( 361 , 361 ) , tt( 361 , 361 ) 
real  phtoh(B,8) 
real  ccos( 59 ) , ssin( 59 ) , c2( 59 ) 

OLD  DATA  FILES 

open(unlt-5,file-’transform.sph' , 6tatus»'old' ) 
open(unit-8,flle-’lcuu.6ph' , status-'old' ) 
open(unlt-9,file-'mas.sph’ , status-’old' ) 
open (unit-10, f ile-'liue.sph' , status- 'old' ) 
open(unit-ll,flle-'l(uel.sph' , status-’old' ) 

NEW  FILE  FOR  STORING  THE  REDUCED  DYNAMICAL  MATRIX 
open ( uni t-4 , f ile-’ reddynmtrx.dat' , status- 'new' , form- 'unformatted' ) 

pi-acos( -1 . 0 ) 

ASSEMBLE  THE  FULL  DYNAMICAL  MATRIX 

Read  in  Ruu(360,360) 
do  i-1,360 

do  j«l, 360,5 

read(  8 ,  *  )  (  kuu(  i  ,  l( ) ,  k-j,j+4) 
end  do 
end  do 

Read  in  Muu(360,360) 
do  i-1,360 

do  j-1,360,5 

read(  9 ,  *  )  (mas(  i  ,  l< ) ,  k»j,j  +  4) 
end  do 
end  do 

Read  Kue(360,l)  into  the  361th  column  of  dyn(361,361) 
do  i-1,360 

read( 10 , * )dyn( i , 361 ) 
end  do 

Read  Rue  transpose  (1,360)  into  the  361th  row  of  dyn(361,361) 
do  j-1,360 

read( 11 , * )dyn( 361 , j ) 
end  do 

dyn( 361 , 361 )-capaci tance  for  zero  displacement  everywhere 
dyn(361,361)-1.8515495E-08 

C  ASSEMBLE  R-w**2M  MATRIX 

wr i te( *,*) 'Enter  frequency  in  Hz;’ 
read( * , * ) f req 
omega-(2.0*pi*f req) 
omega2-omega**2 
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do  i-1,360 
do  j«l,360 

dyn( i , j )-kuu( i , j )-omega2*mas( i , j ) 
end  do 
end  do 

TEST  SYMMETRY  OF  ASSEMBLED  FULL  DYNAMICAL  MATRIX 
nerr-0 
c  do  i-1,361 

c  do  j-i,361 

c  i£  (dyn( j , 1 ) .ne .dyn( i , j  ) )  then 

c  write  (*,*)' Symmetry  error  in  full  dyn.  matrix  at  (i,j)-',i,j 

c  nerr-nerr+1 

c  end  if 

c  end  do 

c  end  do 

c  write  (*,*)  'Number  of  symmetry  errors  in  full  dynamical  matrix  -',nerr 

c  if  (nerr.gt.O)  then 

c  STOP 

c  end  if 

c 

c  REARRANGE  ROWS  SO  ZERO  (INTERNAL)  FORCES  COME  FIRST 

do  j-1.361 
r(l, j)-dyn(l, j) 
r(2, j)«dyn(2, j) 
r ( 241, j )-dyn( 3, j  ) 
r(3. j)-dyn(4, j) 
r(4, j)-dyn(5, j) 
rjs, j)»dyn(6, j) 
r (6, j )-dyn( 7 , j ) 
r(242, j)-dyn(8, j) 
r(243, j)-dyn(9, j) 
r(7, j)-dyn(10, j) 
r(8,j)-dyn(ll, j) 
j)-dyn(12, j) 
cll0,j}~dyn(13,j) 
r(244, j)-dyn(14, j) 
r(245, j)-dyn(15, j) 
r(ll, j)=dyn(16, j) 
r(12, j)-dyn(17, j) 
r(13, j)-dyn(18, j) 
r(14, j)-dyn(19, j) 
r ( 246, j )-dyn( 20, j ) 
r(247, j)-dyn(21, j) 
r( 15, j )-dyn( 22, j ) 
r( 16, j )-dyn(23, j  ) 
r(17, j)-dyn(24, j) 
r(18, j)-dyn(25, j) 
r(248, j)-dyn(26, j) 
r(249, j)-dyn(27, j) 
r(19, j)-dyn(28, j) 
r(20, j)-dyn(29, j) 
r(21, j)-dyni30, j) 
t(22, j)-dyn(31, j) 
r(250,J)-dyn{32, j) 
rj  251 , j )-dyn( 33, j  ) 
r(23, j)-dyn(34, j) 
r(24, j)-dyn(35, j) 
r(25, j)»dyn(36, j) 
r(26, j)-dyn(37, j) 
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r(252, j)-dyn(38, j) 
t(253,l).dyn(39, j) 
r(27, j)-dyn(40, j) 
t(28, j)-dyn(41, j) 
r(29, j)-dyn(42,j) 
r(30, j)-dyn(43, j) 
r(254, j)-dyn(44,j) 
r(2b5,j)-dyn(45, j) 
1(31, J)«dyn(46, j) 
r(32, j)-dyn(47, j) 
r(33, j)-dyn(48, j) 
t(34, j)-dyn(49, j) 
t(256, j)-dyn(50, j) 
r(257, j)-dyn(51, j) 
r(35, j)-dyn(52, j) 
r(36, j)-dyn(53, j) 
t(37, j)-dyn(54, j) 
r(38,j)-dyn(55,j) 
r(258,j)-dyn(56, j) 
r(259,j)-dyn(57, j) 
r(39, j)-dyn(58. j) 
r(40, j)-dyn(59, j) 
c(41, j)-dyn(60, j) 
r(42, j)-dyn(61, j) 
r(260, j)-dyn(62, j) 
r(261, j)-dyn(63, j) 
r(43, j)-dyn(64, j) 
r(44, j)-dyn(65, j) 
r(45,J)-dyn(66, j) 
r(46, j)-dyn(67, j) 
r(262, j)-dyn(68, j) 
r(263, j)-dyn(69, j) 
r(47,  j)-dyn(70, j) 
r(48, j)-dyn(71, ji 
r(49.j)-dyn(72,j) 
r(50.j)-dyn(73,j) 
r(264, j)-dyn(74, j) 
t(265, j)-dyn(75, j) 
r(51, j)-dyn(76, j) 
r(52, j)-dyn(77, j) 
r(53,J)-dyn(78, j) 
r(54. j)-dyn(79, j) 
r(266, j)-dyn(80, j) 
r(267, j)-dyn(81, j) 
r(55, j)-dyn(82,j) 
r(56, j)-dyn(83, j) 
r(57, j)-dyn(84,  j) 
r(58, j)-dyn(85, j) 
r(268, j)-dyn(86,  j) 
r(269, j)-dyn(87, j) 
r(59, j)-dyn(88. j) 
t(60, j)-dyn(89,  j) 
r(61,5>-dyn(90,  j> 
r(62, j)-dyn(91,  j) 
r(270, j)-dyn(92,  j) 
r(271,j)-dyn(93,  j) 
r(63, j(-dyn(94, j» 
r(64, j)-dyn(95,  j) 
r(65, j)-dyn(96,  j) 
r(66, j»-dyn(97,  j) 
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r(272,j)-dyn(98. j) 
r(273, j)-dyn(99, j) 
r(67. j)-dyn(100,j) 
r(68. j)-dyn(101. j) 
r(69. j)-dyn(102, j) 
r(70. j)-dyn(103. j) 
r(274, j)-dyn(104. j) 
r(275, j)-dyn(105,j) 
r(71, j)-dyn<106, j) 
r(72, j)-dyn(107, j) 
t(73, j)-dyn(108. j) 
r(74, j)-dyn(109, j) 
t(276, j)-dyn(110, j) 
r(277,  j)-dyn(lH.  j) 
r(75, j)-dyn(112, j) 
r(76, j)-dyn(113, j) 
r(77,  j)-dyn(U4,  j) 
r(78, j)-dyn(115,j) 
r(278,i)-dyn(116, j) 
r(279,j)-dyn(117, j) 
r(79. j)-dyn(118, j) 
r(eO, j)-dyn(119, j) 
r(81, j)-dyn(120, j) 
r(82, j)-dyn(121, j) 
r(280. j)-dyn(122, j) 
r(281, j)-dyn(123, j) 
r(83, j)-dyn(124, j) 
r(84, j)-dyn(125, j) 
r(85, j i-dyn( 126, j  j 
r(86, j)»dyn(127, j) 
ri282, j)-dyn(128, j) 
r(283, j)-dyn(129, j) 
r(87,j)-dyn(130, j) 
r(88, j)-dyn(131, j) 
r(89, j)-dyn(132,  j) 
t(90, j)-dyn(133, j) 
r(284, j)-dyn(134, j) 
t(285, j)-dyn(135, j) 
r(91, j)-dyn(136, j) 
r(92, j)-dyn(137,i) 
r(93. j)-dyn(138,  j) 
r(94, j)=dyn(139,  j) 
r(286, J)-dyn(140,  j) 
1(287. j)-dyn(141,j) 
r(95,j)-dyn(142,  j) 
r(96, j)-dyn(143,  j) 
r(97, j)-dyn(144,  j) 
r(98,  j)-dyn(14';,  j) 
r(288, j)-dyn(146, j) 
t( 289, j l-dyn( 147  ,  j  ) 
r(99, j)-dyn(148, j) 
r(100,j)-dyn(149, j) 
r(101, j)-dyn(150, j) 
r(102, j)-dyn(151, j) 
r(290, j)-dyn(152, j) 
r(291,J)-dyn(153, j) 
r(103, j)-dyn(154, j) 
r(104, j)-dyn(155, j) 
r(105, J)-dyn(156, J) 
r(106,J)-dyn(157,j) 


r(292, 

j)-dyn( 

158, 

j) 

r(293, 

j)-dyn( 

159, 

j) 

r(107, 

j)-dyn( 

160 

j) 

r(108 

J)-dyn( 

161 

j) 

t(109, 

j)-dyn( 

162, 

j) 

r(110, 

j)-dyn( 

163, 

j) 

r(294 

j)-dyiH 

164 

j) 

r(295 

j)-dyn( 

165 

j) 

j)-dyn( 

166, 

j) 

r(112 

j)-dyn( 

167 

j) 

r(113 

j)-dyn( 

168 

j) 

r(114 

J)-dyn< 

169 

1) 

r(296 

j)-dyn( 

170 

j) 

r(297 

j)-dyn( 

171 

j) 

r(115 

j )-dyn( 

172 

j) 

r(116 

j)-dyn( 

173 

j) 

r(117. 

j)-dyn( 

174 

j) 

r(lie 

j)-dyn( 

175 

j) 

r(298 

j)-dyn( 

176 

j) 

r(299 

j  )-dyn( 

177 

j) 

r(119 

j)-dyn( 

178 

j) 

r(120 

j)-dyn( 

179 

j) 

r(121 

j)-dyn( 

180 

j) 

r(122 

j)-dyn( 

181 

j) 

r(300 

j)-dyn( 

182 

j) 

r(301 

j)-dyn( 

183 

j) 

r(123 

j)-dyn( 

184 

j) 

r(124 

j)-dyn( 

185 

j) 

r(125 

j)-dyn( 

186 

j) 

r(126 

j  )-dyn( 

187 

J) 

r(302 

j)-dyn( 

188 

j) 

r(303 

j  )-dyn( 

189 

j) 

r(127 

j)-dyn( 

190 

j) 

f(128 

j)-dyn( 

191 

j) 

r(129 

j)-dyn( 

192 

j) 

c(130 

j)-dyn( 

193 

j) 

r(304 

j)-dyn( 

194 

j> 

r(30S 

j  )-dyn( 

195 

j) 

r(131 

j )-dyn( 

196 

j) 

r(132 

j)-dyn( 

197 

-j) 

r(133 

j )-dyn( 

198 

j) 

r(134 

j)-dyn( 

199 

j) 

r(306 

j)-dyn( 

200 

-j) 

r(307 

J)-dyn( 

201 

-  j) 

r{135 

jl-dyn( 

202 

j) 

r(136 

j)-dyn( 

203 

j) 

r(137 

j)-dyn( 

204 

j) 

r(138 

j)-dyn( 

205 

-j) 

r(308 

j)-dyn( 

206 

j) 

r(309 

j)-dyn( 

207 

j) 

r(139 

j)-dyn( 

208 

j) 

c(140 

j)-dyn( 

209 

-j) 

r(141 

j)-dyn( 

210 

j) 

r(142 

j)-dyn( 

211 

-j) 

r(310 

j)-dyn( 

212 

-j) 

r(311 

J)-dyn( 

213 

-j) 

r(143 

j)-dyn( 

214 

-j) 

c(144 

j )«dyn( 

215 

-j) 

r(145 

j)-dyn( 

216 

-j) 

t(146 

j)-dyn( 

217 

-j) 
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r(352,j>-dyn(338,j) 
r(353, j)-dyn(339,j) 
r(227, j)-dyn(340, j) 
r(228,j)-dyn(341, j) 
r(229,j)-dyn(342. j) 
r(230, j)-dyn(343, j) 
r(354, j)-dyn(344, j) 
r(355, j)-dyn(345, j) 
r(231, j)-dyn(346, j) 
r(232, j)-dyn(347, j) 
r(233,j)-dyn(348. j) 
r(234,j)-dyn(349,j) 
r(356,j)-dyn(350, j) 
t(357,j)-dyn(351, j) 
r(235, j)-dyn(352,j) 
r(236,j)-dyn(353,j) 
r(237,j)-dyn(354, j) 
r(238,j)-dyn(355, j» 
r(358,j)-dyn(356, j) 
r{359,j)-dyn(357, j) 
r(239,j)-dyn(358,j) 
r(240, j)-dyn(359, j) 
r(360. j)-dyn(360, j) 
r(361,j)-dyn(361, j) 
end  do 

REARRANGE  COLUMNS  SO  INTERNAL  DISPLACEMENTS  CONE  FIRST 

do  1-1,361 

g(l,l)-t(i,l) 

g(i,2)-r(l,2) 

g(i,241)-r«i,3) 

g(i,3)-r(i,4) 

g(i,4)-r(i,5) 

g(l,5)-t(i,6) 

g(i.6)-r(i,7) 

g(i,242)-r(i,8) 

g(i,243)-r(i,9) 

g(i,7)-r(i,10) 

g(i,8)-t(i,ll) 

g(i,9)-r(i,12) 

g(i,10)-r(i,13) 

g(i,244)-r(i,14) 

g(i,245)-t(i,15) 

g(l,ll)-c(i.l6) 

g(l,12)-r(i,17) 

g(i,13)-r(i,18) 

g(i,14)-r(l,19) 

g(i,24e)-c(i,20) 

g(i,247)-r(i,21) 

g(l,15)-r(l,22) 

g(i,16)-r(i,23) 

g(i,17)-r(i,24) 

g(i,18)-r<i,25) 

g(l,248)-r(i,26) 

g(i,249)-r(i,27) 

g(i,19)-r(i,28) 

g(i,20)-r(i,29) 

g(i,21)-r(i,30) 

g(i.22)-r(i,31) 

g(i,250)-r(i,32) 
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g(l,271)-r(l,93) 

g(i,63)-r(l,94) 

g(1.64>-r(i,95) 

g(i.65)-t(i,96) 

g(i,66)-r(l,97) 

g(l,272)-r(l,98) 

g(i,273)-r(i,99) 

g(i,67)-r(l,100) 

g(i,68)-c(i,101) 

g(i,69}-c(1.102) 

g(i,70)-r(i,103) 

gU,274)-r(l,104) 

g(l,275)-r(i,105) 

g(l,71)-r(l,106) 

g(l,72)-r(i,107) 

g(l,73)-r(i,10e) 

g(i,74)-t(i,109) 

g(i.276)-r(l,110) 

g(l,277)-r(i,lll) 

g(i,75)-r{i,112) 

g(l,76)-t(l,113) 

g(l,77)-r(i,114) 

g(l,7e)-r(l,115) 

g(i,278)-r(i,116) 

g(i.279)-r(i,117» 

g(l,79)-r(l,118) 

g(i.80>-r(i.ll9) 

g(i,81)-r(i,120) 

g(i,82)-r(i,121) 

g(i,280)-r(i,122) 

g(l,281)«r(i,123) 

g(l,83)-r(l,124) 

g(i,84)-r(i,125) 

g(i,85)-t{i,126) 

g(l,86)-r(l,127) 

g(l,282)T(i,128) 

g(i,283)-r(i,129) 

g(l,87)-r(i,130) 

g(i,88)-r(l,13U 

g(i,89)-t(i,132) 

g{l,90)-r(i,133) 

g(i.284)-r(i,134) 

g(i,285)-r(i,135) 

g(i,91)-r(i,136) 

q(i,92)-rU.137) 

g(i,93)-r(i,138) 

g(l,94)-r(i,139) 

g(l,286)-r(i,140) 

g(i,287)-r(i,141) 

g(l,95|-r(l,142) 

g(i,96)-r(i,143) 

g(i,97)-r(i,144) 

g(i,98)-t(l,145) 

g(l,288)-r(l,14e) 

g(l,289»-r(l,147) 

g(i,99)-r(i,148) 

g(i.l00)-r(l,149) 

g(l,101)-r(i,150) 

g(l,102)-r(l,151» 

g(i,290)>r(i,152) 
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g(l,291)-r(i,153) 

g(i,103)-r(i,154) 

g(i,104)-r(i,155) 

g(i,105)-r(i,156) 

g(i.l06)-r(i,157) 

g(i,292)-r(l,158) 

g(l,293)-r(i,159) 

g(i,107)-r(i,160) 

g(i,10B)-r(i.l61) 

g(i,109)-r(l,162) 

g(l,110)-r(i,163) 

g(i,294)-r(l,164) 

g(i,295)-r(1.165) 

g(i.lll)-r(i,166) 

g(i,112)-r(i,167) 

g(i,113)-r(i,160) 

g(l,114J-r<l,169) 

g(i,296)-r(i,170) 

g(i,297)-r(i.l71) 

g(i,115)-r(1.172) 

g(l,116)-r(i,173) 

g(i.ll7)-r(i,174) 

g(i,118)-r(i,175) 

g(i,298)-r(i,176) 

g(i,299)-t<i,177) 

g(i,119)-r(i.l78) 

g(i,120)-r(l,179) 

g(l,121)-r(i,180) 

g(i,122)-r(i.l8U 

g(i,300)-r<l,182) 

g(i,301)-r(i,183) 

g(i,123)-r(i,184) 

g(i,124)-r(l,185) 

g(l,125)-r(i,186) 

g(i,126)-t(l,187) 

g(i,302).r(i,188) 

g(i,303)-r(l,189) 

g(l,127)-r(i.l90) 

g(i,128)-t(i.l91) 

g(i,129)-r(i.l92) 

g(i,130)-r(i,193) 

g(i,304)-r(i,194) 

g(i,305)-r(i,195) 

g(i,131)-r(i,196) 

g(i,132)-r(i,197) 

g(i,133)-r(i,198) 

g(i,134)-r(i,199) 

g(i,306)-r(l,200) 

g(l,307)-r(i,201» 

g(i,135)-r<i,202) 

g(i,136)-r(i,203) 

g(i,137)-r(i,204) 

g(i,13S)-r(l,205) 

g(i,308)-r(i,206) 

g(l,309)-r(l,207) 

g(i,139)-r(i,208) 

g(i,140)-r(i,209) 

g(l,l41)-r(i,210) 

g(i,142).r(i,211) 

g(l,310)-r(1.212) 


g(i.l<3 

9(1,144 

9(1.145 

9(1.146 

9(1.312 

9(1.313 

9(1,147 

9(1.148 

9(1,149 

9(1,150 

9(1.314 

9(1.315 

9(1.151 

9(1.152 

9(1,153 

9(1.154 

9(1.316 

9(1.317 

9(1.155 

9(1.156 

9(1.157 

9(1.158 

9(1.318 

9(1,319 

9(1.159 

9(1,160 

9(1.161 

9(1.162 

9(1.320 

9(1.321 

9(1,163 

9(1,164 

9(1,165 

9(1,166 

9(1,322 

9(1.323 

9(1,167 

9(1,168 

9(1.169 

9(1.170 

9(1.324 

9(1,325 

9(1,171 

9(1.172 

9(1.173 

9(1.174 

9(1,326 

9(1.327 

9(1.175 

9(1.176 

9(1.177 

9(1.178 

9(1.328 

9(1,329 

9(1,179 

9(1,180 

9(1.181 

9(1.182 

9(1.330 


-r(l 

213 

-r(i 

214 

-r(l 

215 

-r(l 

216 

-r(l 

217 

-r(l 

218 

-r(l 

219 

-r(l 

220 

-r(l 

,221 

-r(l 

222 

-t(l 

,223 

-r(i 

224 

-r(l 

225 

-r(l 

226 

-r(l 

227 

-r(l 

228 

-r(l 

,229 

-r(l 

,230 

-r(l 

231 

»r(  i 

232 

-r(l 

233 

-r(  i 

,234 

-r(i 

235 

-r(i 

236 

-r(i 

237 

-r(  i 

238 

-r(i 

239 

-r(i 

240 

-r(l 

,241 

-r(l 

242 

-r(l 

243 

-r  ( i 

,244 

-r(l 

,245 

-r(l 

,246 

-r(l 

,247 

-t(l 

,248 

-r(l 

,249 

-r(i 

,250 

-r(l 

,251 

-t(  1 

,252 

-r(i 

,253 

•  r(  i 

,254 

-t{i 

,255 

-r  (  i 

,256 

-r(  1 

,257 

-r(l 

,258 

-r(l 

,259 

-r(  1 

,260 

-t(l 

,261 

-r(l 

,262 

-r(i 

263 

-c(l 

264 

-c(l 

265 

-t(i 

,266 

-r(l 

,267 

-r(  1 

,268 

-r(l 

,269 

-r(l 

,270 

-r(l 

,271 

-r(  1 

,272 
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g(i,331)-r(i,273) 

g(i,183)-r(i,274) 

g(i,ie4)-r(l,275) 

g( i,185)-r(i,276) 

g(i,186)-r(i,277) 

g(l.332)-r(i,278) 

g(i,333)-r(i,279) 

g(l,187)-r(i,280) 

g(i,188)»r(i,281) 

g(i,189)-r(l,282) 

g(i,190)-t(i,283) 

g(i,334)-r(i,284) 

g(l,335>-r(i,285) 

g(i,191)-r(i,286) 

g(l,192)-r(i,287) 

g(l,193)-r(i,288) 

g(l,194)-r(l,289) 

g(i,336)-r(i,290) 

g(i,337)-r(i,291) 

g(i,195)-t(i,292) 

g(l,196)-r(i,293) 

g(i,197)-r(i,294) 

g(i,198)-r(i.295) 

g(i,338)-r(i,296) 

g(i,339)-r(i,297) 

g(l,199)-r(i,298) 

g(i,200)-r(i,299) 

g(i,201)-r(i,300) 

g(i,202)-r(i,301) 

g(i,340)-r(i,302) 

g(l,341)-r(i,303) 

g(i,203)-t(i,304) 

g(l,204)-r(l,305) 

g(l,205)-r(i.306) 

g(i,206)-r(i,307) 

g(i,342)-t(I ,306) 

g(l,343)-c(i,309) 

g(i,207)-r(i,310) 

g(i,208)-r(i,311) 

g(l,209)-r(i,312) 

g(i,210)-r(i,313) 

g(i,344)-r(i,314) 

g(i,345)-r(i,315) 

g(i,211)-r(i,316) 

g(i,212)-r(i,317) 

gU,213)-r(i,318) 

g(i,214)-r(i,319) 

g(i,346)-c(i,320) 

g(i,347)-r(i,321) 

g(i,215)-r(l,322) 

g(i,216)-r(i,323) 

g(l,217)-r(l,324) 

g(i,218)-r(i,325) 

g(i,348)-c(l,326) 

g(i,349)-r(i,327) 

g(i.219)-r(l,328) 

g(i,220)-t(i,329) 

g(i,221)-r(i,330) 

g(i,222)-r<i,331) 

g(i,350)-r(i,332) 
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g(l,351)-r(l,333) 
g(l,223)-r(i,334) 
g(l,224)-r(l,335) 

9(1.225)-r(i.336) 
g(i,226)-t(i,337) 
g(i,352)-t(i,338) 
g(l,353)-r(i,339) 
g(i,227)-r(i,340) 
g(l,228)-r(i,341) 
g(l,229)-r(i,342) 
g(l,230)-r(i.343) 
g(i,354)-t(i,344) 
g(i,355)-r(i,345) 
g(i,231)-r(i,346) 
g(i,232)-r(l,347) 
g(i,233)-r(1.348) 
g(i,234)-r(i,349) 
g(i,356)-r{i,350) 
g(l,357)-r(l,351) 
g(l,235)-r(l,352) 
g(i,236)-r(1.353) 
g(i,237)-r(i,354) 
g(i,238)-r(i,355) 
g(i,358)-r(i,356) 
g(i.359)-r(i,357) 
g(i,239)-r(l,358) 
g(i,240)-r(i,359) 
g(i,360)-r(l,360) 
g(l,361)-r(i,361) 
end  do 

ASSEMBLE  TRANSFORMATION  MATRIX,  t( 361, 361) 

WHICH  CONVERTS  SURFACE  DOFs  TO  NORMAL  AND  TANGENTIAL 

c  Upper  left  corner 

do  i>l,240 
do  j-1,240 

t(l,j)-0.0 
end  do 
t(i,i)-1.0 
end  do 
c 

c  Upper  right  corner 

do  i-1,240 

do  j-241,361 

t(i, j)-0.0 
end  do 
end  do 
c 

c  Lower  left  corner 

do  i-241,361 
do  j-1,240 

t(l, j)-0.0 
end  do 
end  do 
c 

c  Lower  right  corner 

do  i-241,361 
do  j-241,361 
t(i, j)-0.0 
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end  do 
end  do 

t(241, 2411— 1.0 
t(360,360)-1.0 
t(361,361).1.0 

READ  IN  SURFACE  NODE  COORDINATES  AND  COMPUTE  SIN  THETA  AND  COS  THETA 
FOR  ALL  BUT  FIRST  AND  LAST  SURFACE  NODES 
NOTE  THETA  IS  THE  POLAR  ANGLE 

do  1-1,59 

read(5, ’ (2£10.8) ' )  xx.yy 
rr-sqtt( xx**2+yy**2 ) 
ssin(  1  )*=xx/rr 
ccosj i )-yy/rr 
end  do 

k-241 

do  1-242,358,2 

t( 1 , 1 )-ssln( 1-k ) 
t(  1 , 1  +  1 )— ccos(  i-k  ) 
t( 1  +  1 , 1  )-ccos( i-k ) 
t(i+l,i+l)-ssin(i-k) 
k-k+1 
end  do 

ASSEMBLE  tt(  361,361)  -  TRANSPOSE  OF  t( 361, 361) 
do  1-1,361 
do  j-1,361 

tt( j , 1 )-t( 1 , j  ) 
end  do 
end  do 

MULTIPLY  g  BY  t 
do  1-1,361 
do  j-1,361 

suin-0 . 0 
do  k-1,361 

sum-sum+g (i,k)*t(k,j) 
end  do 

pint( 1 , j )-sum 
end  do 
end  do 

do  1-1,361 
do  j-1,361 

sum-O . 0 
do  k-1,361 

sum-sum+tt( i,k)*pint(k,j) 
end  do 

prod( 1 , j )-sum 
end  do 
end  do 

TEST  SYMMETRY  OF  TRANSFORMED  FULL  DYNAMICAL  MATRIX 
nerr-0 
c  do  1-1,361 

c  do  j-i,361 

c  if  (prod( j , 1 ) .ne .prod( 1 , j ) )  then 

c  write  (*,*)'Syiiim  error  in  transf  full  dyn  mtrx  at  (i,j)-',i,j 

c  nerr-nerr+1 
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c  end  if 

c  end  do 

c  end  do 

c  write  (*,*)  'Number  of  symmetry  errors  in  transf  full  dyn  mtrx  «',nerr 

c  if  (nerr.gt.O)  then 

c  STOP 

c  end  i f 

c 

c  REARRANGE  ROWS  SO  ZERO  (TANGENTIAL)  SURFACE  FORCES  CONE  FIRST 

do  1-1,240 
do  j-1,361 

s(i, j)-prod(i,l) 
end  do 
end  do 
do  j-1,361 
s( 300, j ) -prod ( 241 , j ) 
s(301, j)-prod(242, j) 
s( 241 , j )-prod( 243 , j ) 
s(  302, j )«prod( 244 , j  ) 
s(242, j) -prod (245, j) 
s( 303 , j )-prod( 246 , j ) 
s( 243 , j ) -prod ( 247 , j ) 
s(304, j)-prod(248, j) 
s( 244 , j ) -prod ( 249, j ) 
s( 305, j ) -prod ( 250, j ) 
s( 245, j )-prod( 251 , j ) 
s( 306 , j ) -prod ( 252 , j ) 

8( 246 , j ) -prod ( 253 , j  j 
s( 307 , j ) -prod ( 254 , j ) 
s( 247 , j )-prod( 255 , j ) 
s( 308 , j )-prod( 256 , j ) 
s( 248, j )-prod( 257 , j ) 
s( 309, j )-prod(258, j ) 
s(  249 , j i-prod( 259 , j  j 
s(310, j) -prod (260, j) 
s(250, j)-prod(261, j) 
s( 311 , j )-prod( 262, j ) 
s( 251 , j ) -prod ( 263 , j ) 
s(312, j)=ptod(264, j) 
s( 252, j ) -prod j  265, j ) 
s(313, j)-prod(266, j) 
s(253, j)=prod(267, j) 
s(314, j) -prod (268, j) 
s( 254 , j )-prod( 269 , j ) 
s( 315, j )-prod( 270, j ) 
s(255, j) -prod (271, j) 
s(316, j)-ptod(272, j) 
s(256, j)-prod(273, j) 
s(317, j) -prod (274, j) 
s(257, j)-prod(275, j) 
s(318, j)-prod(276, j) 
s( 258, j  j -prod ( 277 , j ) 
s( 319 , j ) -prod ( 278, j  j 
s(259, j)-prod(279, j) 
s(320, j)-prod(280, j) 
s(260, j)-prod(281, j) 
s(321, j)-prod(282, j) 
s(261, j)-prod(283, j) 
s(322, j)-prod(284, j) 
s(262, j)-prod(285, j) 


47 


8(323, J)-prod(2B6,j) 
8(263, j)-prod(287,j) 
8(324, j)-prod(288,j) 
8(264, j)-prod(289,j) 
8(325, j)-ptod(290,j) 
8(265, j)-prod(291,j) 
8(326, j)-prod(292,j) 
8(266, j)-prod(293,j) 
8(327, j)-prod( 294, j) 
8(267, j)-prod(295,j) 
8(328, j)-pcod(296,j) 
8(268, j)-prod(297,j) 
8(329, j)-prod(298,j) 
8(269, j)-prod(299,j) 
8(330, j)-pcod(300,J) 
8(270, j)-ptod(301,j) 
8(331, l)-ptod(302,j) 
8(271, j)-ptod(303,j) 
8(332, j)-prod(304,j) 
8(272, j)-prod(305,j) 
8(333, j)-prod(306,j) 
6(  273, j )"prod( 307, j  ) 
8(334, j)-ptod(308,j) 
8(274, j)-pcod(309,j) 
8(335, j)«ptod(310,j) 
s(275, j ) -prod ( 311 , j ) 
8(336, j)»prod(312,j) 
8(276, j)-prod(313,j) 
8(337, j)-prod(314,j) 
8(277, j)-prod(315,j) 
s(338,l)-prod(316, j) 
6(278, j)-ptod(317,j) 
8(339, jl-prod(318,jl 
8(279, j)-prod(319,j) 
6(340, j)-ptod(320,j) 
6(280, j)-prod(321,j) 
6(341, J) -prod (322, j) 
s(281  -prod(323,j) 

s(342,3.-ptod(324, j) 
8(282, j)-prod(325,j) 
8(343, j)-prod(326,j) 
6(283, j)-prod(327,j) 
8(344, j)-prod(328,j) 
6(284, j)-prod(329,jl 
8(345, j)-prod(330,j) 
6(285, j)-ptod(331,j) 
8(346, j)-prod(332,j) 
8(286, j)-ptod(333,j) 
8(347, j)-prod(334,j) 
6(287, j)-prod(335,j) 
6(348, j)-prod(336,j) 
6(288, j)-prod(337,j) 
8(349, j)-ptod(338,j) 
6(289, j)-prod(339,j) 
8(350, J)-ptod(340,j) 
6(290, j)-prod(341,j) 
6(351, j)-ptod(342,j) 
8(291, J)-prod(343,jl 
6(352, j)-prod(344,j) 
6(292, j)-prod(345,j) 
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8(353, j)-prod(346,j) 

8(293, j)-prod(347,J) 

8(354, j)-prod(348,j) 

8(294, j)-prod(349,j) 

8(355, J)-prod(350,j) 

8(295, j)-prod(351,j) 

8(356, j)-prod(352,j) 

8(296, j)-pcod(3S3,j) 

8(357, j)-ptod(354,j) 

8(297, j)-prod(355,j) 

8(358, j)-prod(356,j) 

8(298, j)-prod(357,j) 

8(359, j)-prod(358,j) 

8(299, j)-prod(359,j) 

8(360, j)-pcod(360,j) 

8(361, j)-prod(361,j) 
end  do 

c 

c  REARRANGE  COLUMNS  SO  TANGENTIAL  SURFACE  DISPLACEMENTS  COME  FIRST 

do  i-1,361 
do  j-1,240 

mat(i, j)»s(l, j) 
end  do 
end  do 
do  i-1,361 
mat(i,300)-s(l,241) 
mat(i,301)-s(i,242) 
fflat(l,241)-s(i,243) 
niat(i,302)-s(i,244) 
mat(l,242)-s(l,245) 
fflat(i,303)-s(i,246) 
mat(i,243)-s(i,247) 
niat(i,304)-s(l,248) 
mat(i,244)-s(i,249) 
mat(i,305)-8(i,250) 
mat(i,245)-s(i,251) 
fflat(i,306)-s(l,252J 
mat(i,246)-s(i,253) 
inat(i,307)-s(i,254) 
inat(i,247)-s(i,255) 
mat(i,306)-s(i,256) 
mat(l,248)-8(i,257) 
mat(i,309)-s(i,256) 
mat(i,249)-s(i,259) 
mat(i,310)-s(i,260) 
mat(i,250)-s(i,261) 
mat(i,311)-s(i,262) 
mat(i,251)-s(i,263) 
nat(i,312)-s(i,264) 
mat(i,252)-8(l,265) 

Mat(l,313)-s(i,266) 
mat(i,253)-s( i,267) 
iiiat(l,314)-s(  1,268) 
iiiat(i,254)-s(  1,269) 
mat(i,315)-8(i,270) 
mat(i,255)-s(i,271) 

«at(i,316)-s(i,272) 

mat(l,256)-8(i,273) 

inat(l,317)-s(l,274) 

inat(l,257)-s(i,275) 
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*at(i.318)-s(i,276) 

Riat(i,2S8)>s(i,277) 

l■at(i,319)-s(i,278) 

*iat(i,259»-s(i.279> 

iBat(i,320)«s(l,280) 

iiiat(i,260)-s(l,281) 

iRat(i,321)-s(l,282) 

•iat(i,261)-s(i.283) 

iiiat(l,322)-s(l,284) 

niat(l,262)-s(i,285) 

niat(i,323)>6(i,286) 

mat(l,263)-s(i,287) 

inat(i,324Us(i,288) 

iHat<i.264)-s(i,289) 

inat(i,325)-s(i,290) 

niat(l,265)-s(i,291) 

i»at(i,326)-e(i,292) 

niat(1.266)-s(l,293) 

«at(i,327)-s(i,294) 

«at(i,2e7)-s(i,295) 

i«at(l,328)-s(i.296) 

inat(i,268)-s(l,297) 

mat(i,329)-6(i,298) 

«iat(i,269)-s(l,299) 

pat(l,330)-s(i,300) 

Riat(t,270)-s(i,301) 

«iat(i,331)-s(i.302) 

i«at(i,271)-6(i,303) 

iiiat(i,332)-s(i,304) 

maid, 272)-s(i, 305) 

ffat(l,333)-s(i,306) 

maid, 273)-s(i, 307) 

iaatd,334)-sd,308) 

i«atd,274)-s(i,309) 

iiiatd,335)-sd.310) 

«atd,275)-sd.311) 

inatd,336)-sd,312) 

matd,276)-sd,313) 

ffatd,337)-sd,314) 

(iiatd,277)-sd,315) 

inatii,338)-s(i,3l6) 

matd,278)-sd,317) 

l•atd,339)-8d,318) 

matd,279)-6d,319) 

matd,340)-sd,320) 

mat (1,280) -8 (1,321) 

mat(l,341)>8(l,322) 

nat(l,281)-8(l,323) 

i«at(l,342)-8(l,324) 

inat(l,282)-8(l,325) 

i«at(l,343)-8(l,326) 

matd,283)-s(l,327) 

iiiat(l,344)-8(l,328) 

Hiat(l,284)-8(1,329) 

mat(l,345)-8(l,330) 

inat(l,285)-8(l,331) 

nat(l,346)-8(l,332) 

i»atd,286)-8(l,333) 

iiiat(l,347)-8(l,334  ) 

iiat(l,287)-8(l,335) 
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«at(i,34e)-s(i,336) 

aiat(i,288)-s(l,337) 

*at(i,349)-s(i,338) 

•iat(i,289)-s(l,339) 

l■at(l,350)-6(i.340) 

mat(i,290)-s(i,341) 

mat(i,351)-s(i,342) 

iiat(l,291)-8(i,343) 

■iat(l,352)-B(1.344) 

iaat(l,292)-s(i,345) 

■iat(i,353)-s(l,346) 

iiiat(i,293)-s(i,347) 

mat(i,354)-s(l,348) 

■at(l,294)-s(l,349) 

■iat(i,355)-s(i,350) 
aiat(i,  295  )-s(  1,351) 
iiiat(i,356)-s(i.352) 
mat(i,296)-s(i,353) 
nat(i,357)-s(1.354) 

Mat(i,297)-6(i,355) 
inat(i,358)-8(  1,356) 
mat(l,298)-s(l,357) 
nat(l,359)>8(i,358) 
mat(l,299)>8(i,359) 
mat(l,360)-s(l,360) 
mat(l,361)>8(l,361) 
end  do 
c 

c  PARTITION  mat  MATRIX  INTO  a,b,c  AND  d  MATRICES 

do  1-1,299 
do  j-1,299 

a(i,  j)-i«at(i,  j) 
end  do 
end  do 
c 

do  1-1,299 

do  j-300,361 

b(l, j-299)-mat(l, j) 
end  do 
end  do 
c 

do  1-300,361 
do  j-1,299 

c(l-299. j)-mat(l, j) 
end  do 
end  do 
c 

do  1-300,361 
do  J-300,361 

d(l-299, j-299)-roat(l, j) 
end  do 
end  do 

COMPUTE  INVERSE  OF  MATRIX  a 
n-299 
np»299 
do  1-1, n 
do  j-l,n 

yd.  j)-0.0 
end  do 
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y(i.l)-1.0 
end  do 

call  ludcinp(a,n,np,indx,d) 
do  j"l,n 

call  lubksb( a,n,np, lndx,y( 1 , j ) ) 
end  do 

MULTIPLY  C  BY  y 
do  i-1,62 
do  j-1.299 

sum-0.0 
do  k-1,299 

sum-sum+c(i,k)*y(k, j) 
end  do 
p(i,  j)-sum 
end  do 
end  do 

MULTIPLY  p  BY  b 
do  1-1.62 
do  j-1,62 

sum-O.O 
do  k-1,299 

sum-sum4p( 1 , k ) *b( k , j ) 
end  do 
pp(l. j)-sum 
end  do 
end  do 

COMPUTE  d-pp 
do  1-1,62 
do  j-1,62 

nn(l, j)-d(l, j)-pp(l, j) 
end  do 
end  do 

TEST  SYMMETRY  OF  REDUCED  DYNAMICAL  MATRIX 
nerr-O 
do  1-1,62 
do  j-1,62 

f racer r-Abs(nn( j , 1 )-nn( 1 , j ) )/Sqrt(AbB(nn( j , 1 )*nn( 1 , j ) ) ) 

If  ( fracerr .ge. 0.001 )  then 

write  (*,*)' Symmetry  error  in  reduced  dyn.  mtrx  at  (i,j)-',i,j 
nerr-nerr+l 
end  if 
end  do 
end  do 

write  (*,*)  'Number  of  symmetry  errors  in  reduced  dyn.  matrix  »',nerr 
if  (nerr.qt.O)  then 
STOP 
end  if 

THIS  COMPLETES  THE  REDUCTION  TO  NODAL  SURF  VEL  DEG  OF  FREEDOM 

STORE  REDUCED  DYNAMICAL  MATRIX  IN  UNFORMATTED  BINARY  FORM 
do  i-1,62 

write(4)(nn(i, j), j'l>02) 
end  do 
C 

end 
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c 


c 


c 


END  OF  MAIN  PROGRAM 

subroutine  ludcinp(a,n,np,indx,d) 
paraneter(nmax-299, tlny-1 .Oe-32) 
diBenslon  a< 299, 299),  lndx<299),  vv(299) 
d-1.0 
do  i»l,n 
aamax-0 . 0 
do  j-l,n 

i£(abs(a( 1 , j ) ) .9t .aamax )  aanax>abs(a( 1, j ) ) 
end  do 

If (aanax.eq.0.0)  pause  'singular  matrix.' 
vv( i )-l . O/aamax 
end  do 

do  j“l,n 

do  i-l,j-l 

sum>a( i , j ) 
do  k»i,i-i 

sum-sum-a (i,k)*a(k,j) 
end  do 
a( i , j  )-sum 
end  do 
aamax-0. 0 

do  i-j,n 

sum-a( i , j  ) 
do  k-i,j-i 

sum-sum-a (i,k)*a(k,j) 
end  do 
a( i , j  )-sum 
dum-vv( i )*abs(sum) 
if(dum.ge. aamax)  then 
imax-i 
aamax-dum 
end  i  f 
end  do 

i f ( j .ne . imax) then 
do  k-l,n 

dum-a( imax,k ) 
a(imax,k)-a( j,k) 
a( j,k)-dum 
end  do 
d--d 

vv( lmax)-vv( j ) 
end  if 

indx( j )-imax 

if(a( j, j) -eq.O. )  a(j,j)-tiny 
if(j.ne.n)  then 

dum-1 .0/a( j , j ) 
do  i-j+l,n 

a(i,ji-a(i,j) *dum 
end  do 
end  if 
end  do 
return 
end 
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subroutine  lubksb(a,n,np, indx.hi ) 
dimension  a(299,299) ,indx(299) ,hi(299) 

ii-0 

do  i»l,n 

m-indx( i ) 

8um>hi(m) 


hi(m)-hi(i) 
if(ii.ne.O)  then 
do  j-ii»i-i 

sum-sum-a( i , J ) *hi { j ) 


end  do 

else  i£(sum.ne.O. )  then 
ii-i 
end  if 
hi( i )-sum 
end  do 


do  i"n,l,-l 
sum>hi ( i ) 
i£<i.lt.n)  then 
do  j-i+l,n 

sum>sun-a ( i , j ) *hi ( j ) 
end  do 
end  if 

hi( i )-sum/a( i , i ) 
end  do 
return 
end 
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APPENDIX  C 


c  PROGRAM  SPHNTRX:  r^^ads  in  reduced  nodal  dynaaical  Matrix  and  produces 

c  dynaMical  Matrix  in  terns  of  first  seven  Le9endre  polynoMials 

real  pi, sum 
real  xx,yy,rr,ccos(61) 

real  Dnod( 62,62 ) ,P( 61 , 7 ) ,Cinv( 7 , 7 ) ,Pt ( 7,61 ) ,D8ph( 8,8 ) 

OLD  FILES 

open (unit-8, f ile-' transform. sph' , status- 'old' ) 

open (unit-9,  file-' reddynmtrx.dat' .status- 'old' , form- 'unformatted' ) 

NEW  FILES 

open(unit-12, f ile-'Osph.dat' , status-'new' ,form-'unfornatted' ) 
open (unit-1 3 , f lie-' Dsphl . lis' , status-'new ) 
open(unit-14,file-'Osph2.1is' .status-'new' ) 

pi-acos(-l.O) 

READ  IN  COORDINATES  OF  ALL  EXCEPT  FIRST  AND  LAST  SURFACE  NODES 
AND  COMPUTE  COS  OF  EACH  ANGLE. 

Note  first  node  is  at  theta-pi  and  last  node  is  at  theta-0  where 
theta  ia  the  polar  angle 
do  i-1,59 

read(8, ' (2fl0.8) ' )  xx.yy 
rr«sqrt( xx**24yy**2 ) 
ccos( i+l )-yy/rr 
end  do 
ccos( 1  )--l .  0 
ccosj  61  )-l . 0 

UNFORMATTED  READ  IN  REDUCED  DYNAMICAL  MATRIX  Dnod(62,62) 
do  1-1,62 

read (9 ) (Dnod( i , j ) , j»l ,62 ) 
end  do 

COMPUTE  NODAL  DISPLACEMENT  TRANSFORMATION  MATRIX,  P(61,7),  USING  THE 
FIRST  7  LEGENDRE  POLYNOMIALS  EVALUATED  AT  COS(NODE3>  THRU  COS(NODBl63) 
Compute  zeroth  through  sixth  order  Legendre  polynomials 
do  i-1,61 

c2-ccos( i ) **2 
P(l,l)-1.0 
P( 1 , 2 )-ccos( i ) 

P(i,3)-(3.0*c2-1.0)/2.0 
P(i,4)-((5.0*c2-3.0)*ccos(i) )/2.0 
P(i,5)-((35.0*c2-30.0)*c2+3,0)/8.0 
P( i,6)-( ( (63.0*c2-70.0)*c2+15.0)*ccos(i) )/8,0 
P(1 ,7)-( ( (231.0*c2-315.0)*c2+105.0)*c2-5.0)/16.0 
end  do 

COMPUTE  TRANSPOSE  OF  NODAL  DISPLACEMENT  TRANSFORMATION  MATRIX  Pt(7,61) 
do  i-1,7 

do  j-1,61 

Pt(i, j)-P( J.l) 
end  do 
end  do 

COMPUTE  INVERSE  OF  LEGENDRE  POLYNOMIAL  NORMALIZATION  CONSTANTS  Cinv(7,7) 
do  i-1,7 
do  j-1 , 7 

Cinv( i , j )«0 . 0 
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end  do 

Cinv(i,i)-2.0/(2*<i-l)+l) 
end  do 


COMPOTE  UPPER  LEFT  7x7  BLOCK  OF  Dsph( 8 , 8 >-Clnv( 7 , 7 ) *Pt ( 7 , 61 ) *Dnod( 61 , 61 ) 

do  i-1,7 
do  j-1,7 
suni'0. 0 
do  k-1,7 

do  IB-1, 61 
do  n«l,61 

8um-suiB+Clnv(l,k)*Pt(k,m)*Dnod(iB,n)»P(n,  J) 

end  do 
end  do 
end  do 

Dsph( i , j )-sum 
end  do 
end  do 


COMPUTE  UPPER  RIGHT  7x1  BLOCK  OF  Dsph( 8 ,8 )-Cinv( 7 , 7 ) *Pt ( 7 , 61 ) ‘KnodC 61 , 1 ) 
do  i»l,7 
sum-O. 0 
do  j-1,7 

do  k-1,61 

Bum-sum+Cinv(i, j)*Pt( j,k)*Dnod(k,62) 

end  do 
end  do 

D8ph(i,8)-sura 

end  do 


COMPUTE  LOWER  LEFT  1x7  BLOCK  OF 


do  j-1,7 
8um-0 . 0 
do  k-1,61 

suiB-8UJn+Dnod(  62 ,  k )  *P(  k ,  j  ) 
end  do 

Dsph(8, j)-sum 
end  do 


Dsph(e,e)-Knodt(l,61)*P(61,7) 


BLOCKED  CAPACITANCE  IS  UNCHANGED 
D8ph(8,8)-Dnod(62,62) 

STORE  Dsph(8,8} 
do  i-1,8 

write(12)(D8ph(i, j), j-l>8) 
wtlte(13,10)(Dsph(l, j) , j-1,4) 
write(14,10)(Dsph(i, j) , j-5,8) 

end  do 

forma t( 3x, 4el7 .7,//) 
end 
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